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EDQNM closure: “A homogeneous 
simulation to support it”, “A 
quasi-homogeneous simulation to disprove it” 

By J.P. Bertoglio K. Squires 1 2 , and J.H. Ferziger 2 


It is known that two-point closures are useful tools for understanding and predicting 
turbulence. Among the various closures, the Eddy Damped Quasi-Normal Marko- 
vian (EDQNM) approach is one of the simplest and, at the same time, most useful. 
Nowadays, direct numerical simulations (DNS) can provide information that can be 
used to test the validity of two-point theories. It is the purpose of the present work 
to use DNS to validate, or improve upon, EDQNM. 

In the first part of this work (“a homogeneous simulation to support it”) we se- 
lected a case for which EDQNM is known to give satisfactory results: homogeneous 
isotropic turbulence. We then evaluated quantities which may be used to test the 
assumptions of two-point closure approximations: spectral Lagrangian time scales. 
Our goal is to make a careful and refined study to validate (and possibly improve) 
the EDQNM theory. 

The aim of the second part of the work ( “a quasi homogeneous simulation to dis- 
prove it”) is, on the contrary, a test of EDQNM in a much more difficult situation. 
Our purpose is to build a reference case for which EDQNM is likely to give poor 
results. We present an attempt to generate a quasi-homogeneous turbulent field 
containing “organized” structures, by artificially injecting them in the initial con- 
ditions. The results of direct simulations using such initial conditions are expected 
to provide a challenge for EDQNM since this kind of field is simple enough to al- 
low comparisons with two-point theories, but at the same time contains “coherent” 
structures which cannot be expected to be accurately accounted for by closures 
based on expansions about Gaussianity. 


1. Lagrangian spectral times in isotropic turbulence 

1.1 The Assumptions of EDQNM 

In the EDQNM theory the growth of the third order moments is limited by the 
introduction of a damping term in the rate equation for the triple correlations. The 
damping coefficient p(k,p,q) is essentially an inverse time scale for the decay of the 
triple correlation among the three wave-vectors of a triad. It is generally assumed 
that p(k,p,q) is the sum of the inverse time scales for the individual wavenumbers: 

p{k,p,q) + v{k 2 +p 2 + q 2 ) = q(k) + 7/(p) + q{q) 
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and i](k) is specified phenomenologically. 

For flows at Reynolds numbers sufficiently high to contain an inertial subrange, 
the damping coefficient must have the form (Orszag 1970) 

r)(k) = aeh 5 


The form proposed by Pouquet et al(1975) for an arbitrary energy spectrum: 



p 2 E(p)dp ]* + vk 2 


is more commonly used. The value of constant A is generally taken to be 0.355. 

According to this expression the time scale at wavenumber k depends only on 
wavenumbers smaller than k (i.e. only the large eddies influence the damping). 
There is however no reason, on physical grounds, to assume that the effect of the 
smaller scales can be neglected. Indeed, in order to devise a two-point closure com- 
patible with the RNG approach of Yakliot and Orszag, Kraichnan (1987) recently 
introduced a model in which it is assumed that only the small scales affect the 
damping. In fact, neither type of expression can be completely correct, and dur- 
ing his stay at the Center for Turbulence Research (CTR) 1987 summer program 
Kraichnan suggested the following possibility of building a time scale that depends 
on all scales: 


rj(k) = A.[ I* p 2 E(p)dp]* +(u+ v>(k))k' 

Jo 


where i />(&) is given by : 


v>(k) = At 



E{p}dp 
r/(fc) + rj(p) 


where A a and Aft are constants. It can be shown that these two new constants must 
satisfy the following relation: 


A a = A — 2(1— log2)^ 

in order that the new form reduces to the old one in the case of an inertial range. 

During the CTR Summer Program, we tested the new form for the damping by 
introducing it into the EDQNM computation code written by Orlandi. We found 
no significant differences from the results obtained using the classical form, as long 
as A a /A remains larger than 0.5. We noticed a small effect on the skewness, (see 
figure 1); the other quantities remain unaffected. 
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FIGURE 1. Evolution of the skewness in time. Comparison of EDQNM results 
obtained with different formulations for the damping, (a) damping given by Pouquet 
et al (1975); (b) damping suggested by Kraichnan (with A 0 /A = .5). 


1.2 Evaluation of i](k) using direct numerical simulation. 

In order to test the expression for rf(k) used in EDQNM against actual turbulence 
time scales, we have to deduce a spectral time scale from the simulation results. 
Comparisons between EDQNM and the Direct Interaction Approximation suggest 
evaluating the time scale using two-time correlations of the velocity field. The 
proper time scale must be invariant under an arbitrary Galilean transformation 
since the triple correlations must not be affected by convective effects. The time 
scale must therefore be a Lagrangian time scale derived from Lagrangian two-time 
correlations, as suggested by LHDIA. For practical reasons, we used the two-time 
correlations introduced by Kaneda (1981): 

Riji^t'lUx'.t'lt') = (ttf(x,<'|<)uj(x',t'|<')) 

where «j(x, f'|/) is the velocity at time t of the fluid particle whose trajectory passes 
through x at time t 1 . One could also use the two-time correlations defined by 
Kraichnan in the LHDIA or ALHDIA theories, or work with the response tensor, 
but the evaluation of Kaneda’s correlations is simpler. 

Spectral information is obtained by Fourier transforming with respect to the 
initial position of the particles. 

To deduce the Lagrangian correlations from the results of a DNS, one has to 
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follow particle trajectories. This is done using Squires’ code, which was developed 
for the study of particle dispersion. 

1.3 Results 

Two direct numerical simulations of homogeneous isotropic turbulence were used. 
Both used Rogallo’s code. In the first case, a 64x64x64 grid was used. In the second, 
a 128x128x128 grid was introduced in order to increase the Reynolds number. 

We present here only a limited number of results. For a more detailed study see 
the paper by Lee, Squires, Bertoglio, and Ferziger (1987). A study preliminary to 
this work was performed in Lyon using Large Eddy Simulation on a 16x16x16 grid. 

In the DNS, it would have been desirable to place particles at each point of the 
computational grid, but doing so would have increased the cost of the computation 
prohibitively. The number of markers is therefore much smaller than the number 
of grid points. Marker placements were of two types. First, to ensure accurate 
large scale statistics, markers were placed on a 16x16x16 equally spaced grid. This 
provides a set of widely spaced particles. To obtain data on the small scales, a 
second set, consisting of 8x8 lines of 64 particles was used. 

In figure 2a, the Lagrangian correlations obtained from the 128x128x128 run 
are given as functions of the dimensional separation time t t . As expected, the 
correlations corresponding to large wavenumbers decrease rapidly. 

The correlations can also be plotted as functions of the time normalized by the 
candidate time scale; an accurate time scale should collapse the results to a single 
curve, figure 2b shows that, when the classical time scale is used, the collapse is 
good at high wavenumbers but not as good at low wavenumbers. When Kraich- 
nan’s expression is used (with A a = A/2), the collapse is slightly improved at low 
wavenumbers, and remains acceptable at high wavenumbers (figure 2c). 

1.^ Problems and future orientations 

The major problem encountered is the lack of sufficient sample at high wavenum- 
bers. The number of particles used (8x8x64) was probably not large enough to allow 
conclusive results concerning the highest wavenumbers in the simulation. 

In the comparison between Lagrangian and Eulerian correlations, another prob- 
lem appears. As expected, for short times and at high k, the Eulerian correlation 
decreases faster than the Lagrangian one; however, for long times, the trend appears 
to be reversed. In fact , for large values of the separation time, the significance of 
Fourier transforming with respect to the initial positions of the markers must be 
questioned. One could try to use another definition for the spectral time scale in this 
case. However, we believe that information deduced from the small time behavior 
is sufficient for the purposes of our study. 

It would be interesting to extend the present study to anisotropic flows, as almost 
no information is available about the effects of mean velocity gradients or anisotropy 
on the damping. The easiest anisotropic case to investigate would be turbulence 
submitted to uniform solid body rotation since, in this case, steady coordinates (in 
the rotating frame) can be used. Due to the lack of time, we could not investigate 
this case during the summer program. It should be the subject of future work. 
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FIGURE 2. Two-time Lagrangian correlations as a function of the time separa- 
tion at different values of the wavenumber, (a) un-normalized time separation, (b) 
time separation normalized using Pouquet’s formula, (c) time separation normalized 
using Kraiclinan’s formula. 
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In the case of isotropic fields, interest in studying the damping term is quite lim- 
ited, since it is known that the existing formulations lead to satisfactory predictions 
of the decay of isotropic turbulence. This is, to a certain extent, true, (and is the 
reason why our goal is the study of anisotropic turbulence). However we believe 
that, even though limited to isotropic turbulence, the results of the present study 
are valuable. In the study of homogeneous isotropic turbulence containing a gap or 
a peak in the spectrum, for example, the accuracy of the expression for the damping 
is of greater importance. Since the choice of the damping term has an important 
effect on the “localness” of the triadic interaction in k space, a better expression 
for 7 / could provide interesting information concerning localized interactions and, 
possibly, insight about the effect of the presence of coherent structures on the en- 
ergy cascade (see Kraichnan, 1987). Furthermore, we believe that, in large-eddy 
simulation, the influence of the subgrid terms on the Lagrangian time scales should 
be investigated. 

2. A quasi-homogeneous field containing “organized” structures 

2.2 Aim of the study and initial conditions. 

The first part of the work was devoted to the study of a case in which EDQNM 
is known to perform well: isotropic turbulence. The aim of the second part is 
significantly different. It consists of an attempt to simulate a case which, although 
simple enough to be handled by closures with reasonable computation times, is 
likely to provide a severe test of EDQNM. 

Our goal is to get insight into the ability of two-point closures to account for 
situations in which coherent structures are present. How are the predictions of 
EDQNM affected by the presence of organized eddies? 

The existence of coherent structures in wall bounded flows is well known. Wall 
bounded flows would, however, require a large amount of computational effort to be 
predicted by EDQNM, and at least at first, simpler situations have to be studied. 

The basic idea of the present work is therefore to artificially inject coherent struc- 
tures into a “quasi-homogeneous” turbulent field (quasi-homogeneous meaning pe- 
riodic rather than homogeneous). 

During the summer program we tried to define initial conditions for a simulation 
of a flow containing “coherent” structures in a periodic cubic box. The results of 
such a simulation could provide reference data to test EDQNM. 

In an attempt to create such a turbulent field, an experiment has been devised in 
Lyon. This experiment consists of generating organized structures in grid turbulence 
by using a grid equipped with small rotating propellers. One hundred counter- 
rotating propellers are used in this experiment (one at each node of the grid); see 
figure 3. The results (see Michard et al, 1986) display interesting behavior, in 
particular, a spike appears in the spectrum; see figure 4. 

Large Eddy Simulations on 16x16x16 grids, run in Lyon, using initial conditions 
containing four vortices embedded in a random field, show that some experimental 
tendencies can be predicted, and, to a certain extent, understood (for example, 
the existence of a periodic distribution of anisotropy). However, due to the lack 
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FIGURE 3. Grid equipped with propellers (Michard et al,1986) 



of resolution of this coarse grid, it has not been possible to inject a peak in the 
spectrum in the flow direction. 

The work done at CTR consists in building an initial field that contains such a 
peak. We used a 64x64x64 mesh. A few simulation runs were started to test the 
initial field. After several attempts, it was decided to define the initial conditions 
in terms of the vorticity field. In the 64x64x64 box, four vortex-containing cells, 
corresponding to the wakes of four propellers, were defined. In each cell three 
vortices are introduced, centered on helical lines corresponding to the wakes of the 
three blades of each propeller. In order to ensure zero circulation in each cell, 
another vortex was introduced on the centerline. A representation of this initial 
field is given in figure 5. 

This field was embedded in a random fluctuating field. For this random field, we 
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FIGURE 5. Representation of the initial “coherent structures” used in the simu- 
lation. These are the idealized wakes of the counter-rotating propellers shown in 
figure 3. 



FIGURE 6. Initial one dimensional spectra (vortex field + isotropic turbulence, 
averaged over the entire x\ — 3*3 space). The mean flow is in the X 2 direction. 

Eu, E 2 2 i — E3 3 . 


used the one created by Lee & Reynolds (1985) with N = 64, and Re — 51.13. 
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The initial spectrum, in the k 2 direction, is given in figure 6. A small peak 
appears at the wavenumber corresponding to the blade passing frequency. 

2.2 Results and directions for future work. 

The peak in figure 6 is relatively weak, although the energy contained in the 
vorticity field is not small compared to the total energy of the field. However, it 
must be pointed out that the plotted spectrum is an average over the X,Z plane, 
and one should investigate the behavior of the one-dimensional spectra obtained at 
a given point in the X,Z plane. 

It is worth noting that the peak appears at the blade passing frequency (due to 
the construction of the initial conditions) whereas the experimental spike is found at 
the rotation frequency of the propellers. During the simulation, the computed peak 
evolved only slightly and no subharmonic spike appeared at the rotation frequency. 
The simulations were, however, limited to one hundred time steps due to time 
constraints; longer runs are necessary to see the expected behavior. 

Another experimental feature is that there is an axisymmetric contraction down- 
stream of the grid. This axisymmetric contraction was introduced to amplify the 
intensity of the peak. Without it, the peak would die immediately; in fact, the am- 
plification was much stronger than expected. The axisymmetric contraction should 
probably be included in the simulation. It is hoped that this work will be completed 
in the near future. It could provide a challenge for EDQNM, as preliminary studies 
show that the experimented decay of the peak is overpredicted by the closure. The 
existence of a reference field deduced from the simulation would be highly valu- 
able since precise analysis of the mechanism occurring in the experiment is hard to 
deduce from the measured quantities. 
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